Section: Research Program

Neuroscience & Neuroendocrinology: Regulation of the Gonadotrope axis

Participants : Benjamin Aymard, Frédérique Clément, Mathieu Desroches, Soledad Fernández García, Albert Granados Corsellas, Elif Köksal, Maciej Krupa, Lucile Megret, Sixtine Passot, Marie Postel, Jonathan Touboul, Alexandre Vidal.

This work was mostly undertaken in the framework of the REGATE (REgulation of the GonAdoTropE axis) Inria Large Scale Initiative Action, that focuses on mathematical neuroendocrinology issues applied to the hypothalamo-pituitary-gonadal (HPG) axis.

Controlled conservation laws for structured cell populations

We have studied the theoretical and numerical questions raised by our multiscale model of follicle selection. This is needed to fully exploit the model potential in terms of biological interpretation and to enable us to forecast the ovulation rate according to the different physiological and endocrine scenarios that we have elaborated [40] . To describe the terminal stages of follicular development on a cell kinetics basis and account for the selection process operated amongst follicles, we have previously developed a multiscale model describing the cell density in each follicle, that can be roughly considered as a ($N\times$ 2D) system of weakly coupled transport equations with controlled velocities and source term [10] , [11] . Even if, in some sense, this model belongs to the class of renewal equations for structured populations, it owns a number of specificities that render its theoretical and numerical analysis particularly challenging: weak nonlinearity due to the moment-based formulation of velocities and source term, discontinuities in the (cell-phase dependent) velocities and densities (due to the mitosis event), 2D effects (e.g. shear or waterproof interface). On the theoretical ground, we have obtained rigorous results on the existence and uniqueness of weak solutions with bounded initial data [56] , so that the well-posedness of the model in its most generic formulation is now well established. In the framework of hybrid optimal control, we have proved that there exists an optimal bang-bang control with only one switching time for the optimal ovulatory trajectory, in the case when the density is approximated by Dirac masses [38] , which in some sense generalizes former results obtained in a low-dimensional ODE case [89] . We can also conjecture that every optimal control is a bang-bang control with only one switching time for our PDE case, but the formal proof of it remains to be stated. From the rigorous reduction (exponential convergence in one of the structuring variable) and averaging of the renewal (mitosis) term, we have obtained a simpler system of coupled nonlinear ODEs (corresponding to the zero and first-order moments of the initial PDEs), from which the dynamics of one given follicle can be studied with respect to the pressure exerted collectively by all other growing follicles, in a dynamical game-like framework. On the numerical ground, we have conceived a new method to deal with the discontinuous coefficients [35] and designed a finite-volume scheme implemented on a parallel architecture [84] to overcome some computational difficulties and perform intensive simulation campaigns.

We have also investigated the physiological balance (as well as pathological or genetically-encoded unbalance) between the oocyte growth and proliferation of follicular cells in the earliest stages of follicular morphodynamics, when the very low number of follicular cells excludes the use of a deterministic formalism. To remain in a dynamical framework consistent (in the limit) with PDE renewal equations, we have adapted a stochastic and discrete formalism initially developed in the framework of ecological modeling (e.g. [88] ) to design a stochastic model of early follicular development with its own specificities [39] : 2D population structuring according both to a space variable (distance from the surface of the oocyte) and an age variable (progression along the division cycle), non-zero sized individuals with possible local overcrowding, multiscale formulation (with three interacting scales intricately merged on the dynamical ground).

Dynamical systems and neuroendocrinology

We have previously proposed in [5] , and further analyzed in [4] , a mathematical model accounting for the alternating pulse and surge pattern of GnRH secretion. The model is based on the coupling between two dynamical systems running on different time scales. The faster system corresponds to the average activity of GnRH secreting neurons, which is forced by the slower system that corresponds to the average activity of regulatory neurons. The analysis of the slow-fast dynamics exhibited within and between both systems allows one to explain the different patterns (slow oscillations, fast oscillations and periodic surge) of GnRH secretion both qualitatively and quantitatively.

In an endocrinology-oriented study, we have explained how the dynamics-based constraints imposed on the model parameters amount to embedding time- and dose-dependent steroid control within the model [23] . We then investigated the plasticity of the model and performed in silico experiments inspired from available experimental protocols: luteal deficiency affecting the surge amplitude, surge blockade induced by administration of luteal levels of progesterone during the follicular phase, short-term effects of either progesterone or estradiol bolus administration on the pulse properties.

On the dynamical ground, further exploration of the model has revealed other possible secretion regimes. In particular, during the transition from a surge back to the subsequent pulsatile phase, a pause consisting of small oscillations superimposed on a long-duration pulse may occur. A detailed examination of the pause has revealed that it is shaped by mixed-mode oscillations (MMO) ; the small oscillations are related to the passage of the slow nullcline of the secreting system through a fold point of its fast nullcline. We have computed families of orbit segment undergoing very brutal transitions upon parameter variation in the vicinity of the fold, by appplying pseudo-arclength continuation algorithms (as implemented in AUTO) to one-parameter families of well-posed two-point boundary value problems. We have derived a variety of reductions that allowed us to obtain results both on the local dynamics near the fold (rigorous characterization of small canards and sectors of rotation) and the global dynamics (existence of an attracting unique limit cycle, which is underlain by a return map) [16] .

We have also started to investigate the question of GnRH neuron synchronization on a mesoscopic scale. We have studied how synchronized events in calcium dynamics can arise from the average electric activity of individual neurons, from seminal experiments of calcium imaging performed on embryonic GnRH neurons [116] . Our model reproduces the occurrence of synchronized calcium peaks, superimposed on asynchronous, yet oscillatory individual background dynamics, as well as additional experimental observations (partial recruitment, doublets of synchronization) [50] . Using phase-plane analysis, we have constrained the model behavior so that it meets not only qualitative but also quantitative specifications derived from the experiments, including the precise control of the frequency of the synchronization episodes.

On a data-oriented ground, we have designed an algorithm (DynPeak) for the monitoring of LH (luteinizing hormone) pulse frequency (that mirrors GnRH pulse frequency in many -but not all- cases), basing ourselves both on the available endocrinological knowledge (pulse shape and duration with respect to the frequency regime) and synthetic LH data generated by a simple model [25] (Joint work with Claire Médigue (hormonal data analysis) and Serge Steer (software development)). We have performed the algorithm on different sets of both synthetic and experimental LH time series. We have further diagnosed possible sources of outliers in the series of IPIs which is the main output of the algorithm.

Innovative computational and theoretical tools for slow-fast dynamics

We have extended the study of the recently discovered torus canard phenomenon [98] , that underlies the transition between the spiking and bursting regimes in neuronal models, and can be roughly considered as the combination of a canard phenomenon with a fast rotating dynamics. We have generalized the previous results to a larger class of bursters (such as the classical Hindmarsh-Rose and Wilson-Cowan models), whose bursting regime ends by a slow passage through a fold bifurcation of limit cycles and we have analyzed the underlying bifurcation structure by means of continuation tools [87] , [92] .

We have developed new approaches to compute one-parameter families of isolas, which are isolated bifurcation branches encountered in multiple timescale dynamics, especially in a neuroscience context (e.g. isolas of spiking, bursting or MMO solutions). The main difficultly consists in computing at once an entire isola and continuing it as a single object in the parameter space, despite its inherent instability. We have proposed a new strategy, implemented as a series of Matlab routines [83] , that enables one to perform multiple parallel continuation runs, subject to specific constraints between the different solution branches. Starting from a known (typically stable) solution obtained by direct simulation, our continuation approach combines the discretization of isolas into (possibly numerous) nodes with the use of periodic boundary conditions and a “phase-like” condition generalizing that implemented for the continuation of limit cycles. In addition, the stability of nodes is checked and possible bifurcations undergone by the nodes or isolas are detected in the course of the continuation.

We have investigated the slow-fast behavior of families of limit cycles in piecewise-linear systems approximations of multiple timescale systems, which are known to reproduce the rich dynamical repertoire of smooth systems while being amenable to more direct analysis. We have revisited previous work from the 1990s in order to complete the definition of a “canard cycle” in this context. We have shown that, even in the partial extension (where the fast nullcline is formed by 3 pieces instead of 4 for the entire extension), key features of canard cycles, such as the explosive behavior in parameter space and the shape with respect to the fast nullcline, are preserved [43] .

We have extended our previous work [93] on epsilon-free methods, whose main advantage lies in not assuming the presence of a small parameter. In the case of planar slow-fast systems, the main idea is to associate strong changes of curvature with loci of inflection points of the flow in the phase plane projection, in order to detect transitions from fast to slow epochs and vice-versa and to estimate the timescale ratio when it is hidden. We have shown that inflection lines, that can be easily computed, provide a good approximation to the excitability threshold [7] . We have also studied the possible topological configurations of inflection lines, both in the singular limit and away from it, both in the “canard regime” (where the canard point corresponds to a tangency between two connected components of the inflection set) and in the “relaxation regime”.